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ABSTRACT 


A  spatial  marching  analysis  Is  given  for  economical  computation  of  three- 
dimensional  viscous  subsonic  flows  in  rotating  geometries.  The  governing 
equations  are  based  on  a  small  scalar  potential  approximation  for  the 
vector-decomposed  secondary  flow  velocity.  No  approximation  is  needed  for  the 
streamwise  pressure  gradient  term  and  this  allows  strong  viscous  secondary 
flows,  coordinate  curvature  and  system  rotation  effects  to  influence  these 
pressure  gradients.  This  approach  is  applied  to  three-dimensional  laminar  and 
turbulent  flows  in  rotating  90  degree  hends  and  in  rotating  straight  pipes  and 
ducts.  The  predicted  structure  of  these  flows  is  consistent  with  experimental 
observations  and  measurements.  Computer  solutions  obtained  using  500,000  grid 
points  require  only  about  15  minutes  of  CRAY-1 S  run  time.  This  approach  appears 
promising  for  further  development  and  application  to  centrifugal  impeller  and 
other  turbomachinerv  flows. 


INTRODUCTION 


The  influence  of  system  rotation  on  turbulent  internal  flows  is  known  to 
have  a  major  effect  on  turbomachinery  performance,  vet  the  phenomenon  is  poorlv 
understood.  Although  used  fairly  extensively,  centrifugal  machines  have  in  the 
past  been  designed  largely  on  an  empirical  basis.  In  support  of  the  empirical 
approach,  the  acquisition  of  experimental  data  from  actual  machines  has  been 
quite  extensive.  However,  few  fundamental  experiments  designed  to  isolate 
various  rotation  effects  have  been  performed  and  consequently  the  understanding 
of  the  fundamental  fluid  mechanics  of  centrifugal  devices  is  still  rather 
sparse.  The  objective  of  the  present  investigation  is  to  develop  an  economical 
method  for  computing  three-dimensional  viscous  flows  with  svstem  rotation,  whose 
use  will  enhance  the  understanding  and  allow  prediction  of  important  physical 
effects  of  rotating  flow  in  centrifugal  turbomachines. 

The  extreme  complexity  of  most  turbulent  rotating  flows  presents  an  obstacle 
for  the  utilization  of  computational  methods.  Such  flows  are  three-dimensional 
and  character i zed  hv  large  secondary  vorticitv  and  velocity  generated  by  turnin’ 
and  Coriolis  effects,  turbulent  shear  layers,  tip  clearance  effects  and  other 
multiple  length  scale  flow  structures.  Solution  of  tilt1  three-dimensional 


s-'',  .e-V---, 


averaged  Navier-Stokes  equations  avoids  making  physical  approximations  other 


than  those  associated  with  turbulence  modeling.  However,  this  approach  is  very 


costly,  even  with  modern  supercomputers,  because  the  flow  structures  of 


practical  interest  are  very  complicated  and  accurate  resolution  of  these  flow 
structures  is  expected  to  require  very  large  mesh  densities.  To  avoid  this  high 
cost  of  solution,  physical/mathematical  approximations  have  been  developed  which 
reduce  the  steady  subsonic  Navier-Stokes  equation  to  a  non-elliptic  form  which 
is  well-posed  for  solution  as  a  spatial  forward-marching  initial/boundary-value 
problem.  The  advantage  of  such  an  approach  is  that  forward-marching  solution 
algorithms  can  be  devised  which  are  much  less  costly  in  terms  of  computer 
resources  (run  time  and  storage)  than  algorithms  for  the  elliptic  Navier-Stokes 
equations.  The  trade-offs  are  that  (a)  additional  error  due  to  the  physical/ 
mathematical  approximations  are  introduced,  and  (b)  the  range  of  flow  problems 


which  can  realistically  be  addressed  is  restricted  relative  to  the  Navier-Stokes 


equations  because  of  factors  such  as  flow  separation,  stagnation  points  and 


transonic  effects.  Nevertheless,  this  approach  seems  well  suited  for  a  number 


of  rotating  internal  ducted  flows  at  or  near  design  conditions,  and  can  provide 


very  high  resolution  of  three-dimensional  viscous  flow  structures  at  relatively 


low  cost.  In  addition,  the  spatial  marching  approach  can  provide  a  large  number 
of  detailed  flow  calculations  at  moderate  cost,  for  use  in  design  optimization 


studies . 


Two  basic  types  of  physical/mathematical  approximations  have  been  suggested 


to  reduce  the  Navier-Stokes  equations  to  a  non-elliptic  form  well-posed  for 


forward-marching  solution.  Approximations  in  both  viscous  and  inviscid  terms  in 
the  Navier-Stokes  equations  are  necessary  to  obtain  non-elliptic  (well-posed) 
approximating  eauations.  First,  a  coordinate  system  for  flow  geometry  being 
considered  must  be  constructed  such  that  streamwise  (marching)  coordinates  can 


be  identified.  The  viscous  approximation  entails  neglecting  terms  representing 


streamwise  diffusion.  In  addition,  two  types  of  inviscid  approximations  have 
been  suggested:  (a)  an  assumed  form  for  the  streamwise  pressure  gradient  term, 


and  (b)  a  small  scalar  potential  approximation  for  the  secondary  flow.  F,ither 


of  these  approximations  produces  non-elliptic  governing  equations. 


The  Inviscid  approximation  which  assumes  a  given  form  for  the  streamwise 
pressure  gradient  term  has  obvious  roots  In  two-dimensional  boundary  layer 


theory,  and  has  been  used  extensively.  Variants  of  this  approach  have  been 
used,  for  example,  bv  Patankar  and  Spalding  [Ref.  1],  Carreto,  Curr  and  Spalding 
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[Ref.  2),  Briley  [Ref.  3],  Chia  and  Sokhey  [Ref.  4],  Kreskovsky,  Briley  and 
McDonald  [Ref.  5]  and  Levy,  Briley  and  McDonald  [Ref.  6]  to  compute  flow  without 
system  rotation  and  by  Howard,  Patankar  and  Bordynuik  [Ref.  7]  and  Majumdar, 
Pratap  and  Spalding  [Ref.  8]  to  compute  flows  with  system  rotation.  Although 
the  approximation  of  streamwise  pressure  gradients  can  provide  reasonable 
accuracy  for  many  problems,  this  type  of  approximation  does  not  allow  the 
streamwise  momentum  equation  to  he  influenced  by  experimentally  observed 
distortions  of  the  static  pressure  field  which  are  induced  by  large  secondary 
flows  associated  with  duct  curvature  and  system  rotation. 

A  second  type  of  inviscid  approximation  (termed  the  small-scalar  potential 
approximation)  has  been  investigated  recently  by  Briley  and  McDonald  (Ref.  9). 
This  approximation  does  not  employ  an  approximation  for  streamwise  pressure 
gradient  terms  and  instead  approximates  convective  terms  in  the  secondary-flow 
momentum  equations  by  neglecting  the  scalar-potential  component  of  a  vector- 
decomposed  secondary-flow  velocity  field  which  corrects  the  transverse  velocity 
vector  from  an  a_  priori  potential  flow  solution.  This  approximation  allows 
strong  viscous  secondary  flows  and  curvature  terms  to  influence  streamwise 
pressure  gradients  in  the  primary  flow  momentum  equation.  It  should  be  noted 
that  the  small  scalar  approximation  is  especially  advantageous  for  rotating 
flows  because  such  flows  are  generally  rotational  even  when  inviscid,  and  this 
precludes  any  convenient  method  of  introducing  an  inviscid  pressure 
approximation,  such  as  imposing  streamwise  pressure  gradients  from  a  potential 
flow. 

In  the  present  report,  a  derivation  of  the  approximating  small  scalar 
potential  equations  is  given  for  a  rotating  coordinate  system.  Computed  results 
for  three-dimensional  viscous  flow  in  simple  confined  rotating  geometries  are 
given.  Some  results  for  laminar  flow  in  a  rotating  90  degree  bend  and  for 
turbulent  flow  in  rotating  straight  ducts,  previously  reported  by  Lin,  Briley 
and  McDonald  (Ref.  10)  are  included  as  part  of  the  present  final  report.  In 
addition,  computed  resuLts  for  turbulent  flow  in  a  rotating  HO  degree  bend  are 
given . 

ANALYSIS 

The  sma L L  scalar-potential  approximation  is  considered  in  detail  for 
non-rotating  coordinate  svstems  in  Ref.  9.  The  present  analysis  considers  the 


derivation  of  the  applicable  approximating  equations  for  a  rotating  coordinate 
system.  The  governing  equations  are  derived  through  approximations  made 
relative  to  a  curvilinear  orthogonal  coordinate  system  fitted  to  and  aligned 
with  the  flow  geometry  under  consideration. 

Equations  governing  the  primary  flow  velocity,  l'p,  and  a  -.cnuidurv  flow 
vorticity,  ft,  normal  to  the  transverse  coordinate  surface-  ire  derived  utilizing 
approximations  which  permit  solution  of  the  equations  as  m  iui'i  if  value 
problem.  Terms  representing  diffusion  normal  to  the  trunsvei-.e  inordinate 
planes  are  neglected.  The  contribution  of  scalar-potential  component s  of 
secondary  velocities  in  convective  terms  of  the  cross  flow  momentum  equations  is 
assumed  small  and  is  neglected.  No  approximation  is  made  for  pressure  gradient 
terms,  and  the  static  pressure  field  is  determined  as  part  of  the  forward¬ 
marching  solution  process.  Secondary-flow  velocities  are  determined  from  scalar 
and  vector  potential  calculations  in  transverse  coordinate  surfaces  once  the 
primary  velocity  and  secondary  vorticity  are  known. 


Compressible  Navier-Stokes  Equations  in  Rotating  Coordinate  System 


The  continuity  and  momentum  equations  for  steady  compressible  flow  relative 
to  rotating  coordinate  system  in  vector  forms  are  given  by 


V/3  u  =  o 


p  N1  =  p  {[(U  •  V)u]  +  p~'  Vp  -  F  +  2(w  X  u)  +  w  x  (oj  X  r  )|  =  0  (2) 

where  p  is  density  and  U  is  velocity  relative  to  rotating  coordinate  system.  p 
denotes  pressure  and  F  denotes  force  due  to  viscous  stress.  w  is  the  system 
rotation  vector,  r  is  the  position  vector.  The  equation  of  state  for  a  perfect 
gas  is  given  hv 

p  =  pRT  O) 


where  R  Is  gas  constant  and  T  Is  temperature.  In  the  present  studv,  the  flows 
are  low  Mach  numher  subsonic  flows  with  neglible  heat  transfer,  so  the 


stagnation  energy  can  he  assumed  to  he  a  known  constant,  E0,  and  the  energy 
equation  can  he  omitted  from  consideration-  For  constant  stagnation  energy,  the 


gas  law  can  he  written  as 


p  =  -T"  _U"2) 


where  y  is  the  specific  heat  ratio. 

In  the  following,  coordinates  x,v,z,  velocity  components  u,v,w  and  the  unit 
vectors  i- ,12*13,  in  the  x-,  y-  and  z-  directions,  respectively,  refer  to  a 
general  rotating  orthogonal  coordinate  system.  The  metric  coefficients  are 
denoted  h^,  ho  and  h^. 

Secondary  Velocity  Decomposition 

The  analysis  is  based  on  a  decomposition  of  the  secondary  velocity  vector 
into  vector  components  derived  from  scalar  and  vector  surface  potentials, 
denoted  4>  and  tp ,  respectively.  The  velocity  vector  is  written  as 

U  =  Tu  +  U,  (5) 

where  u  is  the  primary  velocity  and  Us  denotes  the  secondary  velocity.  The 
secondarv  velocitv  is  written  as 


=  '2 v  +  '3* 


and  is  decomposed  as  follows: 


o,  =  +  -r- 


Here,  7S  is  the  transverse  surface  gradient  operator,  given  for  orthogonal 
coordinates  bv 


n  *  I  d  •>  \  d 

v‘  =  +  '3  t;  it 
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The  transverse  velocity  components  can  he  expressed  as 


v 


w  =  w 


* 


+  w, 


(9) 


using  (7),  (8),  (9),  the  decomposed  transverse  velocity  components  are  given  in 
orthogonal  coordinates  by 


i  At  .  i  i 

(10a) 

h2  dy  h(h3  p 

di 

\  d4>  ii 

dh,  ^ 

h3  <3*  h,h2  p 

dy 

(10b) 

The  decomposition  of  secondary  velocity  v,  w  into  v  ,  v  ,  w  ,  w  ,  introduces 

9  <P  9  V 

two  additional  dependent  variables  and  thus  requires  two  additional  equations  to 
close  the  set  of  governing  equations.  The  additional  equations  are  obtained 
from  vector  identities  associated  with  the  decomposition,  as  discussed  later. 


Physical  Approximations 


Inviscid  Approximations  for  Convective  Terms 

For  convenience  in  defining  the  present  approximations,  a  parameter,  S,  (to 
be  assigned  a  value  of  0  or  1)  is  introduced  in  the  expressions  for  transverse 
velocity  components  as  follows, 

v  =  v4>  +  V '  *  =  P  Vf  +  (  1  1 } 

W  =  ,  w  =  (12) 

The  parameter,  B,  wi 1 L  be  used  to  define  approximations  in  the  convective  terms. 

I'sing  (II)  and  (12),  the  components  of  the  convective  term  C(n)t(n*V)H  can 
be  expressed  in  orthogonal  coordinates  as  follows: 


( n) 


C,  =  U-Vu  +  u(  v  K)2  +  »K|3)  -  v2K2|  -  w2K3i) 

c2  =  U  VV  -  u(uK|2  -'vK21)-w(«KJ2  -VKa)  (14) 

C3  =  U-  Vw  -  u{  uK,3  -  wKjj)  +  v{WKk  -  V  K23)  (15) 

where  the  quantities  K-j  ^  are  the  geodesic  curvatures  of  the  coordinates, 

defined  hy 


(h'h.r'|r 


(  K») 


i 

and  in  which  Xj.X2.x3  are  interchangeable  with  x.y.z,  respectively.  The  small 
scalar-potential  approximation  is  made  hy  setting  B=0  in  (11)  and  (12)  and 
hence,  in  C2 , ^  3  ?  no  approximation  is  made  in  Cj.  If  8=1,  then  v,w=v,w,  and  the 
above  express  ions  for  C2,C3  revert  to  C2,C 3,  their  exact  forms. 


Viscous  Approximation  Neglecting  Streamwise  Diffusion 
The  viscous  force,  F,  in  (2)  can  be  written  as 


PT  =  -  VX(Mil)  +  (X  +  2^)V(V  U)  (17) 


where  *1  =  V x|i  is  vorticity  and  y  and  X  are  viscosity  coefficients.  For  moderate 
subsonic  Mach  number  V  •  U  is  small,  and  the  last  term  in  (17)  is  neglected. 
All  viscous  terms  which  contain  either  a  derivative  with  respect  to  x  or  those 
containing  v*.  or  w^  are  of  smaller  order  than  remaining  terms  and  are 
neglected.  Neglecting  these  terms  give  the  following  approximations: 


r<>r  the  viscous  force  Fj' 


(h, 


♦  ±  Alit Ak  4.  a 

dz  h*  dz  'h,'  ^  1  2  0  di 


(18) 
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and  for  the  transverse  components  of  viscous  force  F 2 ’ ,  F3': 


(h(h3)  p?z 


d(h ,  fJ.fi,) 
dz 


(19a) 


(h,h2) 


3y 


(19b) 


where  the  streamwise  vorticity  is  given  by 


ft. 


(h2h3) 


a(h5w^,) 

-  <*y 


a(h2v^) 

dz 


(20) 


System  of  Approximating  Equations 


Introducing  the  inviscid  approximation  for  convective  terms  and  the  viscous 
approximation  neglecting  streamwise  diffusion,  and  combining  equations  (2), 

(13),  (14),  (15),  (18),  (19)  and  (20),  the  compressible  Navier-Stokes  equations 
can  be  approximated  as: 


pM|  -  p  j 

c,  ♦  Wff 

-  F,‘  +  2(wXy), 

+  [w  X  (wX?)]J 

►  =  0 

(21) 

PM 2  =  P  \ 

^  *  ‘eV1  £ 

-  F2  +  2(wXv)2 

+  [w  X  (wX7)]2j 

\  -  0 

(22) 

PM 3  =  P  \ 

[s  *  W  £ 

-  Fj  ♦  2(w  Xv)3 

+  [w  X  (w  X 7 )] 

[  -  ° 

(23) 

The  continuity  equation  and  equation  of  state,  i.e.,  Eqs.  (1)  and  (4)  remain 
unapproximated,  and  the  decomposed  secondary  velocities  and  11^  satisfy 


following  relations 


The  ahove  seven  equations  provide  a  system  of  equations  governing  the  five 


velocity  components  u,  v^ ,  w^,  v^,  w^,  pressure  p  and  density  p. 

Approximate  Equations  as  a  V/ell-Poseri  Initial  Value  Problem 

It  was  shown  in  Ref.  9  that  for  nonrotating  coordinates  the  foregoing 
approximations  produce  a  system  of  equations  which  is  well-posed.  In  rotating 
coordinates,  the  equations  differ  only  in  the  appearance  of  centrifugal  and 
Coriolis  terms.  These  terms  do  not  affect  the  well-posedness  of  the  equations, 
and  consequently,  the  system  of  equations  (1),  (4) ,  (21),  (22),  (23),  (24a) 
and  (24b)  is  well-posed  for  solution  as  an  initial /boundary  value  problem. 

Approximate  Equations  Written  in  New  Dependent  Variables 

Although  the  dependent  variables  and  the  approximate  equations  given  above 
are  convenient  for  the  analysis,  they  are  not  convenient  for  numerical  solution 
The  equations  are  reformulated  for  numerical  solution  in  terms  of  the  axial 
velocity,  u,  pressure,  p,  streanwise  vorticitv,  flj,  scalar  and  vector  surface 
potentials,  $  and  ip,  together  with  density,  p.  The  equations  for  fij  and  p  are 
derived  hv  taking  the  divergence  and  curl  of  the  transverse  vector  momentum 

'S  A 

equations  (22)  and  (23).  Let  Ms  =  i  2  +  i  3*^  3  denote  the  vector  transverse 
momentum  equation,  the  equations  governing  streamwise  vorticitv,  flj,  and 
pressure,  p,  are  given  by 


'i  '  (VX  =  0 


(3r») 


=  0 


(  ?b  ) 


Incorporating  the  definition  of  11,4,  and  H^,,  the  continuity  equation  (1) 
becomes 


V-p(  iu  +  Vf  )  — 


0 


( :/) 


and  the  definition  of  Sij  becomes 

n,  =  f,  •  vx(/>-' vxi,*) 


(28) 


The  final  system  of  equations  consists  of  equations  (2S),  (2f>),  (27)  and  (28) 
above,  the  x  momentum  equation  (21),  the  state  equation  (4)  for  the  dependent 
variables  u,  p,  flj  ,  <j>,  t}/  and  p.  These  equations  are  Riven  for  a  general 
rotating  orthogonal  coordinate  system  in  Appendix  A. 


Turbulence  Model 

In  turbulent  flow  calculations,  an  isotropic  eddy-viscosity  formulation  is 
used  for  Reynolds  stresses  as  follows: 


P  ui  uj 


~  M 


(29) 


and  the  effective  turbulent  viscosity, 
u.  The  turbulent  viscosity  is  related 
mixing  length  distribution 


UT*  is  added  to  the  laminar  viscosity, 
to  mean  flow  variables  by  means  of  a 


..  _  ,2/0=  =il/2 

P T  -  P^o  (2e=  e) 


(30) 


where  e  is  the  mean  flow  rate  of  strain  tensor 


1  =  t[(V7)  +  (Vv)T]  (31) 

The  mixing  length,  U0,  is  determined  from  the  empirical  relationship  of 
McDonald  and  Camarata  [Ref.  11]  for  equilibrium  turbulent  boundary  layers  which 
can  he  written 


*0<7) 


0.09Sb  ton  h 


(32) 


where  Is  the  local  boundary  layer  thickness,  <  is  the  von  Karman  constant, 


taken  as  0.43,  v  is  distance  from  the  wall,  and  is  a  sublayer  damping  factor 
defined  hv 

%  =  Pl/2{y*-y  *)/<j  (33) 

where  P  is  the  normal  probability  function,  y+  =  y ( t/ p) 2/ ( u/ p) ,  x  is  local 
shear  stress,  v+  =  23,  and  Oj  =  8.  y  is  taken  as  the  distance  to  the  nearest 
wall. 

It  is  recognized  that  this  treatment  represents  a  major  simplification  of 
the  representation  of  the  turbulent  transport  in  a  rotating  system.  However,  it 
is  also  recognized  that  turbulent  transport  effects  are  most  significant  near 
the  solid  wall,  and  in  this  area  a  length  scale  varying  with  distance  from  the 
wall  provides  a  reasonable  first-order  estimate.  Further  from  the  wall,  this 
length  scale  variation  becomes  inaccurate,  but  here  the  flow  is  essentially 
inviscid  and  the  errors  in  the  length  scale  specification  appear  less 
important.  At  present,  in  addition  to  turbulence  model  considerations,  there 
are  major  issues  in  adequately  representing  the  convective  processes  and  in 
obtaining  accurate  numerical  solutions  of  the  governing  system.  As  a 
consequence,  the  present  simple  turbulence  model  was  considered  a  reasonable 
starting  point  for  the  present  turbulent  work. 

Numerical  Method 

The  governing  equations  are  replaced  by  an  implicit  finite-difference 
approximation.  Three-point  central-difference  formulas  are  used  for  all 
transverse  spatial  derivatives.  An  analytical  coordinate  transformation  devised 
by  Roberts  [Ref.  22]  is  employed  for  each  transverse  coordinate  direction,  as  a 
means  of  introducing  a  nonuniform  grid  to  concentrate  grid  points  in  the  wall 
shear-layer  regions.  Two-point  backward  difference  approximations  are  used  for 
streamwise  derivatives,  although  this  is  not  essential. 

In  all  the  solutions  reported  here,  no-slip  or  symmetry  boundary  conditions 
are  prescribed,  as  appropriate.  No  boundary  condition  is  required  for  density, 
since  it  is  computed  algebraically  from  the  state  equation.  The  no-slip 
condition,  v  =  u  =  0,  must  be  expressed  in  terms  of  4>,  and  flj.  The  normal 
voloritv  component  is  specified  by  prescribing  <[>  =  0  and  the  Neumann  condition 


component  of  the  no-slip  condition  is  written  as 


I,  (Ut  +  Vt4>  +  />“' V  X?(  =  0  °A) 

where  it  denotes  the  unit  tangent  vector,  and  the  finite-difference  forms  of 
(28)  and  (34)  are  combined  to  provide  a  boundary  condition  relating  at  wall 
to  <j>  and  t|>.  A  Neumann  boundary  condition  for  pressure,  p,  at  a  no-slip  wall  is 
obtained  from  the  normal  momentum  equation  as 

i„-[Vp-^F‘+  />wX(w*7)]  =  0  (35) 

where  the  convective  terms  and  Coriolis  force  term  vanish  because  II  =  0.  It  is 
beneficial  to  introduce  a  further  change  of  variables  expressing  the  pressure, 
P,  as 


p  =  P(x,  y,  2)  +  pm(x)  +  Ap(x,y,z) 


(36) 


where  P(x,y,z)  is  the  potential  pressure.  Equation  (26)  now  governs  Ap,  but 
does  not  contain  pra,  which  is  in  effect  the  arbitrary  constant  of  the  Neumann 
problem  for  Ap  at  each  x-location.  As  a  consequence,  equation  (26)  can  be 
solved  for  Ap  before  pm  is  known,  and  pra  can  then  be  adjusted  during 
solution  of  the  x-momentum  equation  to  ensure  that  the  integral  mass-flux 
relation 


constont 


(37) 


is  satisfied,  where  A  denotes  cross  section  area. 

A  summary  of  the  procedure  used  to  advance  the  solution  a  single  streamwise 
step  to  the  (n+l)-level  xn+l  from  known  quantities  at  xn  follows.  Unless 
specifically  mentioned  to  the  contrary,  the  transverse  velocities  v^,  w^, 
v,ij ,  w,p  and  the  density  p  are  evaluated  explicitly  at  the  n-level.  In 
addition,  the  convective  operator  is  evaluated  as  pn!ln  »V.  Values  of  u>,  hj, 
h2  ,  and  hj  are  given  and  thus  known  at  both  xn  and  xI1+'. 

1.  Equations  (A.l)  and  (A. 2)  form  a  linear  coupled  system  for  ft;™  'and  (j/1 
which  is  solved  as  a  2x2  coupled  system.  For  this  purpose  artificial  time 
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derivatives  are  added  to  each  equation,  and  an  iterative  block-implicit  scheme 
[Ref.  12]  is  used.  In  prescribing  no-slip  boundary  conditions,  the  tangential 
component  (34)  contains  a  contribution  from  this  contribution  is  evaluated 

using  4>n.  Terms  in  the  vorticity  equation  (A.l)  containing  u,  v,  w,  and 

nnn\l>y  _  ih  ib 

w  are  evaluated  using  u  ,  v  ,  w  ,  v  and  w  ;  x-deri vati ves  of  v  and  w  are 

evaluated  using  n-  and  (n-l)-level  quantities,  w2,  w3  are  known  and 

evaluated  at  n+1  level. 

(2)  The  pressure  equation  (A. 3)  is  solved  for  Apn+1  using  an  iterative 
scalar  ADI  scheme.  In  this  equation  all  appearances  of  v^+*  and  Wy+^  are 
evaluated  using  ^n+l,  and  u,  v,  w,  p,  v^  and  w^  are  evaluated  using 
n-level  quantities.  r^  ,  r2,  r3  and  ojj,  w2,  to3  are  known  and  evaluated  at  n+1 
leve 1 . 

(3a)  Using  an  assumed  value  of  p^  ^  to  begin  a  secant  iteration  and  values 
of  Apn  and  Apn+1  now  available,  the  x-momentum  equation  (A. 4)  is  solved  to 
determine  un+^ ,  using  a  scalar  ADI  scheme. 

(3b)  The  density  pn+*  is  evaluated  from  the  state  equation  (A. 5)  using 

-n+1  n+1  ,  n+1  ,  .  ,  . .  ,, 

U  ,  Pn,  and  Ap  ,  which  are  now  available. 


(3c)  For  internal  flows  the  integral  mass-flux  relation  (37)  is  evaluated 
using  un+1  and  pn+l. 

(3d)  Assuming  that  the  initial  guess  for  p^+1  was  not  exact,  the  integral 
mass-flux  relation  will  not  be  satisfied,  and  steps  (3a-c)  are  repeated 
iteratively  using  the  standard  secant  method  [Ref.  13]  to  find  the  value  of 
pj^+^  that  leads  to  un+^  and  pn+^  satisfying  the  integral  mass-flux 
relation  (37). 

(4)  Finally,  the  continuity  equation  (A.b)  is  solved  for  4>n+l  using  an 

iterative  scalar  ADI  scheme  and  currently  available  values  of  un+^  and  pn+l. 

,  .  n+1  n+1  n+1  n+1  ,  ,  c  n+1 

The  velocity  components  v  ,  w  ,  v  and  w  are  then  evaluated  from  <j> 

and  ’in+l  . 


COMPUTED  RESULTS  FOR  ROTATINO  FLOWS 

One  of  the  motivations  of  the  present  studv  is  the  eventual  computation  of 
flow  in  centrifugal  turbomachinery  components  such  as  centrifugal  impellers. 

The  treatment  of  radial  impeller  geometries  requires  that  the  present  analysis 
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he  generalized  for  use  with  nonorthogonal  coordinate  systems.  This  generali¬ 
zation  and  the  implementation  In  impeller  geometries  Is  part  of  an  ongoing 
investigation  being  performed  under  a  related  contract  (DAA029-85-C-0030).  In 
the  present  study,  the  analysis  is  applied  to  rotating  flows  in  simple 
geometries  using  orthogonal  coordinates.  The  flows  considered  include  laminar 
and  turbulent  flow  in  a  rotating  90  degree  bend  with  square  cross  section  and 
turbulent  flows  in  rotating  straight  pipes  and  ducts. 

Laminar  Flow  in  a  Rotating  90  Degree  Bend 

The  geometry  and  flow  configuration  considered  are  shown  in  Fig.  1.  The 
flow  enters  the  duct  axially  and  leaves  radially.  The  Reynolds  number 
Re  =  U0 D/v  is  790,  where  UQ  is  the  mean  flow  velocity,  D  is  the  duct  width, 
and  v  is  kinematic  viscosity.  The  duct  has  a  square  cross  section,  and  the  Mach 
number  was  taken  as  0.001,  which  means  that  the  flow  considered  is  essentially 
incompressible.  Both  rotating  and  nonrotating  flows  are  considered  for  this 
geometry.  The  initial  boundary  layer  thickness  6^/D  was  taken  to  be  either 
0.4  or  0.1.  In  Ref.  9,  predictions  for  nonrotating  flow  with  6j/D  =  0.4  were 
found  to  be  in  very  good  agreement  with  experimental  measurements.  For  the 
rotating  flow,  two  cases  are  reported  here,  one  with  Rossby  number  Ro  equal  to 
5.0  and  the  other  with  Ro  equal  to  1.0.  The  Rossby  number  is  defined  as 
Ro  =  'J0/2ftD,  where  fi  is  the  rotation  speed. 

The  small-scalar  potential  approximation  is  examined  in  Fig.  2  for  both 
non  rot  a t ing  flows  and  rotating  flows  with  Ro  =  1.  The  initial  boundary  layer 
thickness  was  0.4.  The  maximum  absolute  value  of  each  ( <(>  or  ip)  component  of 
transverse  velocity  is  shown  against  streamwise  distance.  The  results  in  Fig.  2 
show  that  the  ^-component  of  the  transverse  velocity  is  significantly  less  than 
the  ^-component ,  except  very  near  the  entrance  section  of  the  bend.  This 
comparison  is  taken  as  an  indication  that  the  assumption  of  small  11^  is 
reasonable  for  these  flows.  The  small-scalar  potential  approximation  has  proved 
adenuate  in  all  the  test  cases  considered. 

Computed  results  for  both  rotating  and  nonrotating  flow  cases  with 
5 | / n  =  0.4  are  compared  in  Figs.  3-5.  Shown  in  these  figures  are  the 
transverse  seconda rv-f low  velocity  vectors,  the  streamwise  primarv-flow  velocity 
contours,  and  the  static  pressure  contours  for  the  flow  at  the  00  degree 
location,  at  the  end  of  the  bend.  In  the  nonrotating  flow,  the  curved  duct 
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go one try  generates  a  strong  secondary  flow  as  shown  in  Fig.  3,  and  this  distorts 
the  primary  flow  velocity  as  shown  in  Fig.  4.  The  static  pressure  is  shown  in 
Fie.  5.  The  rotating  flow  cases  are  also  shown  for  comparison  in  Figs.  3-5. 

The  Coriolis  forces  in  the  rotating  flow  exert  an  increasing  influence  on  the 
flow  development  as  the  system  rotation  rate  increases.  The  rotation  causes 
additional  secondary  flow  motion  which  tends  to  twist  the  counter-rotating 
vortices  present  in  the  nonrotating  case,  and  this  further  distortion  increase's 
wi,h  increasing  system  rotation  rate,  as  shown  in  Fig.  3.  The  streamwise 
velocity  is  shown  in  Fig.  4,  and  the  observed  distortion  of  the  primary-flow 
velocity  contours  in  the  rotating  flow  case  is  consistent  with  the  predicted 
secondary  flow.  Finally,  the  static  pressure  in  Fig.  5  reflects  a  complicated 
interaction  of  the  primary  and  secondary  flows  as  influenced  bv  duct  curvature, 
centrifugal,  and  Coriolis  force  effects. 

Fig.  6  and  Fig.  7  show  the  effects  of  initial  boundary  layer  thickness  on 
the  development  of  primary  flow  in  both  nonrotating  and  rotating  90  degree 
bends.  The  thicker  houndary  laver  causes  a  stronger  secondary  flow  which  in 
turn  generates  stronger  distortion  in  the  primary  flow  as  shown  in  Fig.  b  and 
Fig.  7.  Phvsicallv,  the  secondary  flow  is  generated  inviscidly  hv  the  turning 
of  transverse  vorticitv  hut  its  growth  is  retarded  by  viscous  stress  near 
walls.  Viscous  stresses  exert  a  stronger  retarding  effect  on  thin  boundary 
laver  cases  than  on  thick  boundary  layer  cases.  Thus,  the  thicker  boundary 
laver  flow  generates  stronger  secondary  flow  and  distortion  in  the  primary 
flow.  The  computed  effect  of  initial  boundary  layer  on  the  development  of 
primary  and  secondary  flow  is  consistent  with  experimentally  observed  trends 
f  3  e  f .  14], 

Turbulent  Flow  in  a  Rotating  90  Degree  Bend 

The  geometry  and  flow  configuration  are  the  same  considered  for  laminar 
flow  (Fig.  1).  The  flow  parameters  are  Re  =  50,000,  Ro  =  1.0  and  initial  shear 
laver  thickness  S ^ / f»  =  0.1.  This  calculation  was  started  in  the  straight 
extension  upstream  of  the  bend  at  a  distance  of  1.5  D  from  the  start  of  the 
He’d.  The  development  of  this  flow  at  streamwise  stations  corresponding  to  0, 
If',  HO  and  9ft  degrees  of  turning  are  shown  in  Figs.  Ra-d.  The  development  nt 
the  primarv  and  secondary  flow  is  shown  in  Figs.  Ra-b.  At  the  start  of  the  bend 
(*•  degrees)  there  is  little  Coriolis  effect  because  the  primarv  flow  ve loci t  v  is 
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parallel  to  the  axis  of  rotation.  There  is  only  a  small  amount  of  secondary 
flow,  and  this  Is  due  to  the  effect  of  geometric  curvature.  As  the  bend  angle 
of  turning  Increases  (30,  60,  90  degrees)  there  Is  an  increased  effect  of 
Coriolis  forces,  combined  with  secondary  flows  generated  by  turning  of  the 
duct.  The  Coriolis  forces  cause  additional  secondary  flow  which  twists  the 
vortices  induced  by  turning.  The  vortices  migrate  toward  the  forward  wall  in 
the  direction  of  rotation  (the  low  pressure  surface)  and  this  behavior  is 
consistent  with  general  experimental  observations.  A  considerable  amount  of 
distortion  of  the  primary  flow  is  present,  as  is  consistent  with  predicted 
secondary  flows.  Because  the  Coriolis  force  depends  on  the  local  flow  speed, 
the  overall  flow  structure  is  quite  complicated.  The  secondary  flow  velocity  is 
verv  large  and  of  the  same  order  of  magnitude  as  the  primary  velocity  in  thin 
shear  layers  along  the  inside  and  outside  walls  of  the  bend.  These  large 
secondary  vortices  are  the  result  of  the  primary  flow  momentum  deficit  in  these 
shear  layers  subjected  to  Coriolis-induced  pressure  gradients. 

The  pressure  coefficient  for  this  flow  is  shown  in  Fig.  Be.  At  the  start 
of  the  hend ,  the  transverse  pressure  gradients  are  relatively  mild.  The 
gradient  is  essentiallv  in  the  direction  from  outside  to  inside  of  the  bend  and 
is  associated  with  geometric  curvature  effects.  As  the  bend  turning  angle 
increases,  the  transverse  pressure  gradients  gain  considerable  strength  and  are 
directed  essentially  in  the  direction  of  rotation  with  low  pressure  on  the 
leading  surface.  Finally,  contours  of  the  streamwise  pressure  gradient 
(3Cp/9x)  are  shown  in  Fig.  8d .  These  gradients  are  of  interest  because  the 
present  analysis  and  governing  equations  were  formulated  to  avoid  approximations 
for  this  term  in  the  primary  flow  momentum  equation.  These  streamwise  pressure 
gradients  gain  considerable  strength  and  complexity  with  increasing  bend  turning 
angle;  the  transverse  variation  in  this  term  is  of  order  unity.  Tt  would  be 
difficult  to  devise  a  pressure  approximation  for  this  term  which  would  be 
adequate  for  rotating  flow  structures  of  this  complexity. 

Turbulent  Flow  in  a  Rotating  Pipe 

The  geometrv  and  flow  configuration  for  this  case  are  shown  in  Fig.  9.  The 
Reynolds  number  Re  =  ll^d / v  Is  6  X  10^,  based  on  pipe  diameter,  d,  and  mean 
axial  velocity,  !Tn.  The  Rossby  number  Ro  =  nm/2ftd  is  1.0,  where  ft  is  the 
rate  of  rotation.  The  initial  boundary  layer  thickness  is  0.01  d,  which  was 


estimated  from  the  entrance  axial  velocity  profile  measured  by  Kikuyama, 
Murakami  and  Nishibori  [Ref.  15]  for  these  flow  conditions.  Kikuyama,  Murakami 
and  Nishibori  [Ref.  15]  measured  time-mean  velocity  and  turbulent  fluctuations 
in  this  flow.  They  found  the  system  rotation  of  the  pipe  has  a  stabilizing 
effect  in  the  suppression  of  turbulence.  Using  the  mixing-length  theory,  they 
suggested  the  following  formulas  for  mixing  length  distribution  to  account  for 
the  suppression  of  the  turbulence  due  to  pipe  rotation. 


=  I  -  BRi 


where  l  and  in  denote  the  values  of  the  mixing  length  for  rotating  and 
nonrotating  flow.  8  is  a  constant,  and  R|  is  the  Richardson  number  defined  'y 


2w  d  f/ du  \z  d  (  w  ' 

’’  =  T*  aT(rwllbr)  +  I'alTT. 


Here,  w  is  the  swirl  velocity  and  u  is  the  streamwise  velocity,  r  and  z  are 
distances  from  the  center  and  the  wall  of  pipes,  respectively.  In  the  present 
calculation,  the  mixing  length  distribution  l0  for  nonrotating  flow  is 
computed  by  using  equations  (32),  (33),  and  the  mixing  length  distribution  for 
turbulent  pipe  flow  is  calculated  by  using  equations  (38),  (39).  The 
coefficient,  8,  is  chosen  to  be  2.  The  turbulent  viscosity  is  obtained  from 
equat i ons  ( 30 ) ,  (31). 

Figures  (10a)  and  (10b)  show  a  comparison  of  computed  and  measured 
streamwise  velocity,  u,  and  swirl  velocity,  w,  components  at  station  x/d  =  28.8, 
The  streamwise  velocity,  u,  is  normalized  by  free  stream  velocity,  Ue ,  and  the 
swirl  velocity,  w,  is  normalized  by  pipe  radius,  ro,  and  rotation  rate,  S3.  The 
distance  from  the  pipe  wall,  rQ-r,  is  normalized  by  the  momentum  thickness, 

9x.  The  computed  and  measured  velocity  distributions  are  In  good  agreement. 


Turbulent  Flow  in  a  Rotating  Duct 


The  geometry  and  flow  configuration  for  this  test  case  are  shown  in 
Fig.  11.  Tlie  rectangular  cross  section  has  a  height  to  width  ratio  ii/D  =  D.388. 
the  Reynolds  number.  Re,  is  bf>,5DO,  based  on  passage  hydraulic  diameter  and  mean 
velocity.  The  Rossbv  number,  Ro,  is  b.  based  on  passage  width,  D  and  mean 


velocity.  The  initial  boundary  layer  thickness  used  in  the  calculation  is 
n.05  D.  The  mixing  length  distribution  from  equation  (32)  was  used  for  these 
calculations.  In  Figs.  (I2a)  and  (12b),  the  computed  streamwise  velocity  in  the 
horizontal  center  plane  and  secondary  velocity  in  the  vertical  center  plane  are 
compared  at  x/D  =  5  with  the  measurements  of  Wagner  and  Velkoff  [Ref.  16]  and 
the  computed  results  of  Howard,  Patankar  and  Bordynuik  [Ref.  7].  The  present 
results  agree  very  well  with  the  measured  streamwise  velocity  in  Fig.  (12a),  and 
the  improved  agreement  relative  to  the  Ref.  7  calculation  probably  reflects  the 
tiner  grid  and  resolution  of  the  viscous  sublayer  region  in  the  present  results, 
compared  with  the  wall-function  treatment  of  Ref.  7.  The  comparison  for 
secondary  velocity  in  Fig.  (12b)  is  less  conclusive  because  of  the  difficulty  of 
measuring  small  velocities  near  walls  using  hot  wire  anemometry. 


CONCLUSIONS 


A  spatial  marching  analysis  is  described  for  computation  of  three- 
dimensional  viscous  flows  with  system  rotation.  The  governing  equations  are 
based  on  a  small  scalar  potential  approximation  for  the  secondary  flow,  which 
does  not  entail  any  approximation  of  the  pressure  gradient  terms.  The  analysis 
was  assessed  bv  application  to  laminar  and  turbulent  rotating  flows  in  simple 
ducts  and  bends.  This  approach  seems  well  suited  for  a  number  of  rotating 
internal  ducted  flows  at  or  near  design  conditions,  and  can  provide  very  high 
resolution  in  three  dimensions  at  relatively  low  costs.  Adequate  resolution  of 
turbulent  flows  including  viscous  sublayer  resolution  appears  to  require  at 
least  500,000  grid  points  even  in  the  simple  90  degree  bend  geometry.  Solutions 
obtained  by  the  present  method  using  500,000  grid  points  require  onlv  about  15 
minutes  of  CRAY-1S  CPU  time.  Future  work  will  address  implementation  of  the 
analysis  in  nonorthogonal  coordinates  and  computation  of  centrifugal  impeller 
f lows . 
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APPENDIX  A 


Governing  Equations  in  Rotating  Orthogonal  Coordinates 
The  governing  equations  for  the  reduced  form  of  the  Navier-Stokes  equations 
based  on  the  small-scalar  potential  approximation  is  given  helow  In  a  general 
rotating  orthogonal  coordinate  system. 

Vorticity  Equation 
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( A.  1  ) 


Here,  Dj  is  streamwise  vorticity,  and  uij,  u)2,  w3  are  components  of  the  system 
rotation-rate  vector.  Kj ^  is  given  in  equation  (16). 

Vorticity  Definition 
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where  C2  and  C3  .ire  given  in  equations  (17),  (18),  and  r  ;,  r2,  r3  are 
components  of  the  position  vector. 


x-Momentum 
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Figure  Be  -  Development  of  Pressure  Coefficient  Cunt  ours  for  Turbulent  Flow 
in  a  Rotating  90°  Bend. 
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